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I. INTRODUCTION 

The nature of phase transitions in frustrated systems 
are generally found to be quite different from those in 
conventional magnets. [Q Frustration often leads to non- 
trivial ground state degeneracies as in the case of the 
antifcrromagnetic spin 1/2 Ising model on the trian- 
gular lattice. The classical Heisenberg model on 
the two dimensional triangular lattice with antiferromag- 
netic nearest neighbour coupling and easy axis exchange 
anisotropy is another example where frustration leads to 
a novel ground state degeneracy. Miyashita and Kawa- 
mura [|}(MK) have investigated both the ground state 
and the nature of harmonic excitations at low tempera- 
tures in this model. At zero temperature the spins lie in 
a plane which contains the easy z-axis. In addition to the 
Si degeneracy related to the rotation of this plane about 
the z-axis, there is also a non-trivial Si degeneracy of 
the ground state related to rotations of the spins within 
the plane. The planar spin configuration distorts as it is 
rotated about an axis normal to the plane but the energy 
remains constant. This degeneracy is not broken by the 
harmonic excitations. Using Monte Carlo methods, MK 
measured the specific heat as well as various susceptibili- 
ties and suggested that the system undergoes two succes- 
sive phase transitions as the temperature is lowered. The 
two transitions indicated the onset of power law correla- 
tions of the spin-spin correlation functions parallel and 
perpendicular to the easy axis. More recently Sheng 
and Henley || have examined both the effects of finite 
temperature and quantum fluctuations on this non-trivial 
degeneracy. They predicted that the continuous degen- 
eracy at finite temperatures T is reduced to a discrete 
six-fold degeneracy and that an additional phase transi- 
tion from a floating phase to a locked phase occurs at a 
temperature in between the upper and lower transitions 
for large values of the easy axis anisotropy. 

In the present work, we use a direct method to study 
the free energy cost of both spin wave and vortex for- 
mation in this system at low temperatures. We define 
a spin stiffness and a vorticity modulus in terms of an 
equilibrium correlation function which can be evaluated 
using standard Monte Carlo methods. Our results indi- 
cate that there are only two transitions and that both 
correspond to vortex unbinding transitions accompanied 
by power law decay of spin correlations. 



II. THE MODEL 



The Hamiltonian of the system is given by 



H = 



i<j 



r jX qX qX i ry QV qV , jZ qZ qZ] 
[J Oj ~v J Oa ~r o >3j\ 



(1) 



where Sf, a = x,y, z represents a classical 3-component 
spin of unit magnitude located at each site j of a trian- 
gular lattice and the exchange interactions J x = J v — 



J, J z — AJ are restricted to nearest neighbour pairs of 
sites. We shall consider the case where the parameter 
range A > 1 represents an easy axis anisotropy. The 
limit A — > 1 corresponds to the isotropic Heisenberg 
model j(| whereas the limit A — > oo corresponds to an 
infinite spin Ising model. Southern and Xu || have 
previously studied the Heisenberg limit and found direct 
evidence for a vortex unbinding transition at finite tem- 
perature even though the spin stiffness vanishes on large 
length scales. This latter fact indicates that the spin cor- 
relations decay exponentially at all finite temperatures 
T. This behaviour is quite different from that which oc- 
curs in the corresponding two-component xy model || 
where power law decay of the correlations appear at low 
temperatures. 

The triangular lattice can be decomposed into three 
sublattices A, B and C as shown in figure 1 with the 
three spins on each triangle labelled as Sa, Sb and Sc- 
A chirality vector jl(| for each upward pointing triangle 
is defined as follows 



K A = 



3%/3 



(S A x Sb + Sb x Sc + Sc x S A ) (2) 




FIG. 1. The decomposition of the triangular lattice into 
three sublattices A, B and C. The sublattice labels on all 
black triangles are the same and represent only one of the six 
possible choices. 

In the present case of the easy axis antiferromagnet, 
the low temperature spin configuration on each triangle 
has the spins lying in a plane which includes the z-axis 
and which is perpendicular to the chirality vector. We in- 
troduce a local spin space coordinate system for each up- 
ward pointing triangle using the three unit vectors Ka , z 
and J_a = z x Ka- Hence the vectors z and ±a lie in the 
spin plane and Ka is the local normal to the plane. The 
z is a global axis for all triangles whereas the Ka are only 
all aligned at T = 0. As T increases, the local normal di- 
rections fluctuate and the spin plane develops curvature. 
Hence local coordinates are needed to properly define the 
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response properties of the system. 

Information about the rigidity of the system against 
fluctuations can be obtained from the spin wave stiffness 
coefficient. The spin stiffness (helicity) tensor is given 
by the second derivative of the free energy Jfl],[l^] with 
respect to the twist angle about a particular direction in 
spin space. We apply a twist about each of the local axes 
defined above and the corresponding average stiffness is 
given by 



Pa = 



— ((^(e 4 , -u)[ 



J' 3 + ,P 



][S?iq-S7Sf]) 2 ) (3) 



where the superscipts take the values a =J_a, z, K/\ and 
a, /3, 7 are to be taken in cyclic order. S" denotes the 
component of the spin at site i in the average direction of 
the corresponding local unit vectors of the upward point- 
ing triangles to which the spin belongs, are unit vec- 
tors along neighbouring bonds and u is the direction of 
the twist in the lattice. 

We employ compact clusters, as shown in figure 2, con- 
taining N spins with periodic boundary conditions ap- 
plied that are consistent with the sublattice structure. 
The number of spins is related to the linear size L of the 
cluster as N = 3L 2 . 




FIG. 2. Compact clusters of linear size L containing 
TV = 3L 2 spins used in the Monte Carlo calculations are in- 
dicated by the solid lines. In each case periodic boundary 
conditions are used which are compatible with the three sub- 
lattice structure. 

We have used a single spin-flip heat bath algorithm |b| 
to update the spin directions at each Monte Carlo step 
and all thermal averages are replaced by time averages. 
For the largest value of the system size L = 60 studied we 
discard the first 1.5 x 10 4 Monte Carlo steps and perform 
averages over the next 5 x 10 5 steps. 

The values of the stiffnesses at zero temperature can 
be evaluated using the exact classical ground state cor- 
relation functions. We find these values to be 



Pz/J = 



2A + 1 
2(1 + A) 2 



P±/J = A 2 



A 



2(1 + A) 2 

Pk — Pz + p±_ 



(4) 



These three stiffnesses satisfy a perpendicular axis the- 
orem which is consistent with the ground state being a 
planar spin arrangement. In the Heisenberg limit we have 
p z = p± but for large values of A they have different 
energy scales, p z /J ~ i whereas p±/J ~ 4- These 
limiting values can also be obtained from the harmonic 
excitation spectrum, p] 




FIG. 3. The three stiffnesses p z , p± and pK in units of J as 
a function of T/J for system sizes L — 10, 20, 30, 40, 60 when 
A = 2. The dashed curve represents the line p = (2/V)T. 

Figure 3 shows the results obtained for the three stiff- 
nesses as a function of T/J for various system sizes L 
with A = 2. As the temperature increases towards 
T C1 = (0.231 ± 0.005) J, the stiffness corresponding to 
the smaller energy scale, p z , decreases abruptly and a 
strong dependence on system size L is evident indicating 
that above this temperature the system has no rigidity 
against twists of the spins about the easy axis. MK sug- 
gested that this transition corresponds to the unbinding 
of vortices associated with the chirality vector which lies 
primarily in the xy plane. We will provide further sup- 
port for this interpretation in section IV. 

The remaining two stiffnesses do not vanish at T Cl but 
do exhibit a rapid decrease followed by a more steady 
decrease until a higher temperature T C2 — (0.75±0.03) J. 
At this higher temperature, the strong dependence of the 
stiffness px on system size L indicates that the system 
loses rigidity against rotations about the local chirality 
axes. The Kosterlitz-Thouless(KT) theory for the 

xy model in two dimensions predicts a universal value for 
the ratio p c /T c equal to 2/V. The dashed line in figure 3 
corresponds to p/T = 2/tt and intersects the stiffnesses at 
temperatures where finite size effects first appear. The 
stiffnesses p z and px indicate two temperatures where 
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rigidity of the system about the z-axis is first lost followed 
by a loss of rigidity about the local chirality axis with an 
apparent universal value for the ratio p c /T c in each case. 
The third stiffness p± has a change in curvature at the 
two transitions but remains finite to high temperatures. 
Similar behaviour is found for other values of A as well. 




FIG. 4. The spin stiffnesses (a) pK and (b) p z in units of 
J as a function of T/J for system sizes L = 10,20,30 with 
A = 2,4, 6, 8. The dashed curve in both cases represents the 
line p = (2/n)T. 

Figures 4(a) and 4(b) show our results for px and p z 
respectively for the values A — 2,4,6,8. The dashed 
line plotted in both figures corresponds to p/T = 2/tt 
and in each case this line intersects the stiffnesses at the 
temperatures where finite size effects first become signifi- 
cant. These results indicate that the ratio p c /T c appears 
to have the universal value 2/tt for both stiffnesses for all 
values of A > 1. We have used this criterion to estimate 
the values of T Cl and T C2 . 

In the same way that the spin stiffness is a measure of 
the response of the spin system to a twist over the length 
of the lattice, a vorticity || can be defined as the response 
of the spin system to an imposed twist about a given 
axis a in spin space along a closed path which encloses a 
vortex core. This is essentially the response of the system 
to an isolated vortex and can be calculated as the second 
derivative of the free energy with respect to the strength 



of the vortex, or winding number m, evaluated at m = 0. 
We obtain the following expression 



i<3 



^^(E(^)[^][^-^]) 2 > (5) 



4ttT 



where is the distance of site i from the vortex core and 
(j)i is tangent to the circular path in the lattice passing 
through the site i and enclosing the vortex. Here a, (3, 7 
are defined in the same way as for the stiffnesses and 
indicate the axis of rotation of the vortex. 

The V a contain both a core contribution and a part 
which is proportional to InL. By comparing systems of 
different lattice sizes L\ and L2 we can extract the vor- 
ticity modulus v a defined as follows 



using 



V a = C a + v a In L 



V a {L2)-V a {L x ) 



HL2IU 



(6) 



(7) 



where the v a are normalized so that they have the same 
zero temperature values as the corresponding stiffnesses. 
Our approach does not require any change in boundary 
conditions and is applied directly to the antiferromag- 
netic model. 




FIG. 5. The vorticity moduli vk and v z in units of J as a 
function of T/J obtained using various pairs of system sizes 
when A = 2. The solid curves represent the average of all 
pairs for each modulus. The dashed curve represents the line 
v = (2/tt)T. 

Figure 5 shows our results for the vorticity moduli 
obtained by comparing systems of different sizes when 
A = 2. The behaviour of the moduli as a function of in- 
creasing temperature is identical to the spin stiffnesses 
shown in figure 3. The vorticity modulus v± is only 
weakly sensitive to the two transitions and is not plotted. 



4 



The vanishing of the vorticity moduli v z and vk at T Cl 
and T C2 respectively indicates that free vortices appear at 
these two transition temperatures. The results are also 
consistent with a universal value v c /T c = 2/tt. Similar 
behaviour is found for the values A = 2,4, 6, 8. 

The results presented so far do not indicate directly 
that power law decay of correlations are present below 
these two transition temperatures. In the next section 
we present results for the spin-spin structure factor which 
are also consistent with KT transitions. 



III. STRUCTURE FACTOR 

We have also studied the spin-spin structure factor for 
various system sizes L 



S aa (<l) = i- T < S? Sf > e iq (r 



(8) 



where a = x, y, z . In particular, we have studied S = 
S xx + S vy and S zz as a function of both q and T. In both 
cases the Fourier component Q = (47r/3,0) exhibits a 
divergence at and below the two transition temperatures 
T Cl and T C2 respectively. Other values of q only exhibit 
a maximum. This Fourier component corresponds to the 
three-sublattice structure associated with the triangles. 
If we assume power law decay of the spin-spin correlation 
functions, then the structure factors should depend on 
the system size L as follows 



S ZZ (Q)~L 2 -^ 



(9) 



where rj z and r/± are the corresponding correlation length 
exponents. 





L=20-10 




— = L=40-10 




^ L=60-10 




FIG. 6. The correlation function exponents 2 — r]± and 
2 — rj z as a function of T/J obtained by comparing various 
system sizes with A = 2. The dashed line corresponds to 
r, = 1/4. 



Figure 6 shows the values of 2 — r\ z and 2 — r)± as a 
function of T / J obtained by comparing systems of dif- 
ferent sizes with A = 2. The dashed line indicates the 
value 77 = 1/4 predicted by the KT theory and intersects 
both exponents at temperatures where a strong depen- 
dence on system size first appears. Similar results were 
obtained for other values of A. The universal value of 
r/ z = 77 j_ — 1/4 is consistent with the universal value of 
Pc/T c — 2/tt and v c /T c — 2jix obtained in section II. 



T/J 2 - 




T, 



FIG. 7. The transition temperatures T C1 /J (lower) and 
c 2 /J (upper) as a function of A determined using the uni- 



versal value of the stiffness p and the universal value of r\ are 
shown by the squares and circles respectively. The value at 
the Heisenberg point (A = 1) is indicated by the diamond 
and was obtained from previous vorticity results. The solid 
curves are obtained by assuming that these two temperatures 
are proportional to the T — values of p z and p± respectively. 

We have used these two criteria to obtain the tempera- 
tures T cl and T c2 as a function of A as shown in figure 7. 
The values obtained from the stiffnesses and the structure 
factors agree with each other fairly well. The error bars 
in T c reflect contributions from the statistical uncertainty 
of the Monte Carlo measurements, as well as estimated 
systematic contributions due to finite size effects. If we 
make the simple assumption that the dependence of T Cl 
and T C2 on A is the same as the zero temperature values 
of p z and p± respectively, then we have 



T C1 = T c 



T c 



4(2A + 1) 
3(1 + A) 2 
4A 2 (A + 2) 
3(1 + A) 2 



(10) 



where the values of T Cl and T C2 in the isotropic limit 
A = 1 are both equal to T c = (.305 ± .005) J. f§ These 
expressions are plotted as the solid curves in figure 7 and 
the agreement with the data points is remarkable. 

The limit A — > 00 of the present model is equivalent 
to the S — > 00 limit of the antiferromagnetic spin S Ising 
model on the same lattice with coupling AJ. The solid 
curve for T C2 approaches the value T^jAJ ~ .407 ± .007 
in this limit. This value is also in excellent agreement 
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with that obtained by Nagai et al. [[l7| for the S 1 = oo 
Ising model. In the next section we give a more detailed 
microscopic description of these phases. 



IV. MICROSCOPIC PROPERTIES 

The three spins S^,Sb and Sc on each triangle can 
be expressed in terms of the following three vectors 

M = S A + S B + S c 
M fl — Sa — -S B — -S c 



M / = ^(S B -S C ) 



(11) 



The vectors M^j and M/ are the real and imaginary parts 
of the complex Fourier component Mq = + zMj 
where Q = (4^/3,0). In addition, the chirality for each 
triangle can be written as K = §M^ x Mj. 

In the ground state the three spins lie in a plane which 
includes the z-axis. We denote the angle that each sub- 
lattice spin makes with the z-axis by 6 ai and define 



ip = (9 A + e B + e c )/2> 



(12) 



where < ip < 2ir. We measure ip relative to the z- 
axis in the counter-clockwise direction. As discussed by 
Sheng and Henley |||, the nontrivial degeneracy of the 
ground state can be parameterized in terms of the angle 
ip. It has been previously pointed out that both the en- 
ergy and the magnitude of M are constant in the ground 
state manifold |||l8|]. However, it is easy to show that 
the modulus of the complex vector Mq is also constant. 
In the Appendix we give some exact results for the com- 
ponents of these vectors in the ground state. The _L, z 
components of these two vectors depend on ip in quite 
different ways. The components of Mo have periodicity 
27r/3 whereas the vectors M# and M/ have period 2ir. 



A=2 J 



\|/=0 \]/=7i/6 \|/=;i/3 



\|/=7t/2 i|/=2rc/3 



v.. 



\1„ 



A»l 



FIG. 8. Ground state spin configurations for values of ip 
in the range < ip < 2tt/3 for A = 2 and A > 1. The 
horizontal and vertical axes are _L and z repsectively. The 
corresponding directions of ip and Mo are independent of the 
value of A and are indicated in the middle row by the solid 
and dashed arrows respectively. 



Figure 8 shows ground state spin configurations of the 
sublattices in the _L, z plane for values of ip in the range 
< i> < 27r/3. The top row corresponds to A = 2 and the 
bottom row to A — > oo. Also shown are the correspond- 
ing directions of ip and Mo which are independent of the 
value of A. Rotations of the angle ip by 27r/3 correspond 
to interchanges of the sublattice spin directions and to a 
complete rotation of Mo by 2ir. Hence vortex-like config- 
urations of Mo simply correspond to switches of the sub- 
lattice magnetizations. For large values of the easy axis 
anisotropy, the spin configurations become more Ising- 
like. In the range < ip < 7r/3, sublattices A and B 
remain locked parallel to the easy axis while sublattice 
C simply follows Mo- In the range 7r/3 < ip < 2tt/3, the 
roles are interchanged with sublattice A following the di- 
rection of Mo. These sublattice changes as a function of 
ip are similar to the domain transition regions described 
by Wannier (|] for the S — 1/2 Ising model at T = 0. In 
this latter case, the transition regions contribute to the 
macroscopic entropy of the system and to the power law 
decay of correlations with wavevector Q. || 

The moduli of the complex numbers Mq = Mjj + iMf 



and Mq = 



are not constant in the ground 



state but the combination 



M, 



Q 



zMq- = 



'(4A + 2)(A + 2) eift 



{l + Af 



(13) 



has constant modulus and a phase angle ^ which has 
a leading term linear in ip but with an additional small 
6-fold modulation. At low temperatures the chirality vec- 
tor lies in the xy plane and it can be described by the 
complex number 



K~ + iK„ 



Kl 



K*e 



(14) 



where the modulus is not constant but exhibits only a 
weak cos(6ip) variation in the ground state. These two 
phase angles and ^ describe the complex order pa- 
rameters associated with the transitions at T Cl and T C2 
respectively. 

Our Monte Carlo procedure permits us to take snap- 
shots of the spin configurations at various temperatures. 
In order to separate the topological defects from the con- 
tinuous deformations of the spin configurations, we first 
raise the temperature to some fixed value and allow the 
system to reach equilibrium. We then rapidly quench the 
system to low temperatures and allow the system to ap- 
proach a nonequilibrium configuration. The topological 
defects are metastable at low temperatures and require 
a much longer time to disappear than the continuous de- 
formations. In figures 9 and 10 we show snapshots of the 
spin configuration at very low temperatures obtained by 
heating the system above T C1 and T C2 respectively and 
then rapidly quenching. In both cases, (a) describes the 
spatial variation of the chirality angle O and (b) describes 
the angle ip. 



G 



(a) 



V 
(b) 



FIG. 9. Snapshots of (a) the chirality angle O and (b) the 
angle ip at very low temperature following a quench from just 
above T C1 for A — 2 and L = 20. The dots indicate the 
positions of the black triangles shown in figure 1 and the lines 
indicate the local direction of either O or iji. 

Figure 9 shows snapshots of the spin configuration for 
a cluster of size L = 20 with A = 2 after quenching from 
a temperature of T = 0.5 J, just above T Cl . Vortices asso- 
ciated with the angle are clearly visible in (a) whereas 
only continuous distortions of the angle ip are visible in 
(b) . These vortices in O only persist at low T if the sys- 
tem is first heated to a T > T Cl and then quenched. They 
are metastable configurations at this lower temperature 
but become stable at T Cl . Our results provide direct ev- 
idence for a vortex unbinding transition at T C1 involving 
the chirality vector. 



e 

(a) 



(b) 



FIG. 10. Snapshots of (a) the chirality angle O and (b) 
the angle ip at very low temperature following a quench from 
above T C2 for A = 2 and L = 20. The dots indicate the 
positions of the black triangles shown in figure 1 and the lines 
indicate the local direction of either O or ip. 



Figure 10 shows similar snapshots when the system is 
first heated to a temperature T — 3 J, well above T C2 , 
and then rapidly quenched. Vortices in both and ip 
are clearly visible in this case. The vortices in ip are only 
visible after quenching from above T C2 . Hence the upper 
transition corresponds to a vortex unbinding transition 
associated with the sublattice ordering vector Mq. 

We have also studied how the sublattice magnetiza- 



tions change as we update the spin configurations. In 
the temperature range T < T C2 there is a continuous 
sublattice switching that occurs. The time scale for the 
switching depends both on T and system size L increas- 
ing with L and 1/T. We have not been able to determine 
a scaling form for this time scale. Similar behaviour has 
been reported previously for other frustrated systems on 
the triangular lattice. The antiferromagnetic spin 1/2 
Ising model with both nearest and next-nearest neigh- 
bour interactions jl9],|o) as well as the antiferromagnetic 
S = oo nearest neighbour model |l||21,^2| exhibit this 
phenomena in the temperature range where a KT phase 
occurs. In our case, the switching is due to the wave-like 
variations in ip which can produce "vortex-like" varia- 
tions in the local magnetization vector Mo. In the range 
of A that we have studied, we have not observed a tran- 
sition from sublattice switching to a partially ordered 
phase where one of the sublattices is locked in the +z 
direction, another in the —z direction and the third per- 
pendicular to the z-axis. Sheng and Henley S predicted 
that such a locking transition might occur at low enough 
temperature. 



V. SUMMARY 

We have presented strong evidence for the occurence of 
two distinct Kosterlitz-Thouless types of defect-mediated 
phase transitions in the Heisenberg antiferromagnet on 
the triangular lattice with easy axis anisotropy. A nu- 
merical approach was used to calculate the rigidity of 
the system against both spin wave deformations and the 
formation of free vortices at low temperatures for sev- 
eral values of the easy axis anisotropy. In each case 
a universal value of the stiffness or vorticity modulus, 
Pc/T c — v c /T c — 2/tt, was found. In addition, the 
spin correlation length exponent has the universal value 
7] = 1/4. These are the same values that are predicted by 
the Kosterlitz-Thouless theory for defect unbinding tran- 
sitions in two dimensions. The two phase angles identi- 
fied in equations (13) and (14) in section IV describe the 
complex order parameters associated with the transitions 
at T C2 and T Cl respectively. 

We find that the spin stiffnesses and the correspond- 
ing vorticity moduli behave identically for the easy axis 
case A > 1. In contrast, at the isotropic Heisenberg 
limit, the spin stiffness vanishes at large length scales 
whereas the vorticity moduli are non-zero at low T and 
vanish abruptly at a finite temperature. || Our previous 
work on the xy model || also indicated that the vorticity 
and stiffness behave identically. This corresponds to the 
A — ► — oo limit where there are again two transitions but 
they are very close in temperature with the upper tran- 
sition corresponding to an Ising-like transition and the 
lower to a Kosterlitz-Thouless transition. Recent work 
by Capriotti et. al. [f23[ has also found similar behaviour 
in the range < A < 1. Hence the Heisenberg point is a 
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multicritical point where four phase transition lines meet. 
For A > 1 there are two KT transition lines whereas as 
for A < 1 there is an Ising and KT line. 
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APPENDIX: 

Miyashita and Kawamura (J] were the first to identify 
the nontrivial degenerac y o f the ground state configu- 
ration. Other groups plJlql have studied the effects of 
quantum fluctuations on this degeneracy. One indicator 
of the degeneracy is that the magnetization vector Mo 
has a constant magnitude in the ground state indepen- 
dent of the value of ip. In particular we find the following 
exact relations 

< = -(4-r) sin (3^) 



Ml 



A 
A 



y A+V 
M .Mo = (^-j) 2 



cos(3^) 



(Al) 



where ip is measured counter-clockwise from the z-axis in 
the _L, z plane. A rotation of Mo by 2tt corresponds to 
rotating tp by 27r/3 which is simply a cyclic permutation 
of the sublattices on the triangle. Hence vortices in M 
can be associated with sublattice switching. 

In addition to Mo, the complex vector Mq also has 
constant modulus, 



Mq.M q = Mtj.Mtj + M/.M/ 
4A 2 + 1CL4 + 4 



(1+Af 



(A2) 



The moduli of the components of Mq as well, as the 
chirality, are not constant in the ground state but depend 
on tp as follows 



\m z q \ 2 



I^qI 2 

IKI 2 



3A 2 



6A A 
— + (- 



1 



(1 + ^4) 
6A + 3 ,A 

s5 + (~ 



A + l 
1 



4 r V2A + A 2 



2 cos 2 (3V») 
2 sin 2 (3-0) 



3V3 L (A + l) 2 
V1 + 2A 
+ (1+A) 2 



A+V 

(2 + A) cos 2 (3?/>) 



(l + 2A)sin 2 (3V')] (A3) 
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